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Abstract 

Dislocation nucleation is essential to our understanding of plastic deformation, ductility and 
mechanical strength of crystalline materials. Molecular dynamics simulation has played an im- 
portant role in uncovering the fundamental mechanisms of dislocation nucleation, but its limited 
time scale remains a significant challenge for studying nucleation at experimentally relevant con- 
ditions. Here we show that dislocation nucleation rates can be accurately predicted over a wide 
range of conditions by determining the activation free energy from umbrella sampling. Our data 
reveal very large activation entropies, which contribute a multiplicative factor of many orders of 
magnitude to the nucleation rate. The activation entropy at constant strain is caused by thermal 
expansion, with negligible contribution from the vibrational entropy. The activation entropy at 
constant stress is significant larger than that at constant strain, as a result of thermal softening. 
The large activation entropies are caused by anharmonic effects, showing the limitations of the 
harmonic approximation widely used for rate estimation in solids. Similar behaviors are expected 
to occur in other nucleation processes in solids. 
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Nucleation plays an important role in a wide range of physical, chemical and biological 
processes^-. In the last two decades, the nucleation of dislocations in crystalline solids 
has attracted significant attention, not only for the reliability of microelectronic devices^, 
but also as a responsible mechanism for incipient plasticity in nano-materials^ - — and nano- 
indentatior*ii~— . However, predicting the nucleation rate as a function of temperature and 
stress from fundamental physics is extremely difficult. Because the critical nucleus can be as 
small as a few lattice spacings, the applicability of continuum theory^ becomes questionable. 
At the same time, the time scale of molecular dynamics (MD) simulations is about ten orders 
of magnitude smaller than the experimental time scale. Hence MD simulations of dislocation 
nucleation are limited to conditions at which the nucleation rate is extremely high^^. 

One way to predict dislocation nucleation rate under common experimental loading 
rates^- is to combine the transition state theory (TST)^^ and the nudged-elastic-band 
(NEB) method^. TST predicts that the nucleation rate per nucleation site in a crystal 
subjected to constant strain 7 can be written as 



r TST = ^ exp 



knT 



(1) 



where F c is the activation free energy, T is temperature, and ks is Boltzmann's con- 
stant. The frequency prefactor is u = ksT/h, where h is Planck's constant. Note that 
F C (T, 7) = Eci^f) — T5 C (7), where E c and S c are the activation energy and activation en- 
tropy, respectively. Here we assume the dependence of E c and S c on T is weak, which is later 
confirmed numerically for T < 400K. For a crystal subjected to constant stress a, F C (T, 7) in 
Eq. ([I]) should be replaced by the activation Gibbs free energy G C (T, a) = H c (a) — TS c (a), 
where H c is the activation enthalpy. Because the NEB method only computes the activa- 
tion energy, the contribution of S c is often ignored in rate estimates in solids. Recently, an 
approximation of S c (a) = H c (a)/T m is used^, where T m is the surface disordering tempera- 
ture. This approximation was questioned by subsequent MD simulations^. The magnitude 
of S c remains unknown because none of the existing methods for computing activation free 
energies^— has been successfully applied to dislocation nucleation. 

We successfully applied the umbrella sampling^ method to compute the activation free 
energy for homogeneous and heterogeneous dislocation nucleation in copper. Based on this 
input, the nucleation rate is predicted using the Becker-Doring theory^. Comparison with 
direct MD simulations at high stress confirms the accuracy of this approach. Both F C (T, 7) 



and G C (T, a) show significant reduction with increasing T, corresponding to large activation 
entropies. For example, 5 C (7 = 0.092) = 9 ks and 5' c (cr = 2GPa) = 48 k,B are observed in 
homogeneous nucleation. We found that SJ^y) is caused by the anharmonic effect of thermal 
expansion, with negligible contribution from the vibrational entropy. The large difference 
in the two activation entropies, AS C = S c (a) — 6^(7), is caused by thermal softening, which 
is another anharmonic effect. Similar behaviors are expected to occur in other nucleation 
processes in solids. 




FIG. 1. (a) Schematics of the simulation cell. The spheres represent atoms enclosed by the critical 
nucleus of a Shockley partial dislocation loop, (b) Shear stress-strain curves of the Cu perfect 
crystal (before dislocation nucleation) at different temperatures. 



For simplicity, we begin with the case of homogeneous dislocation nucleation in the bulk. 
Even though dislocations often nucleate heterogeneously at surfaces or internal interfaces, 
homogeneous nucleation is believed to occur in nano-indentatior*ii and in a model of brittle- 
ductile transition^. It also provides an upper bound to the ideal strength of the crystal. Our 
model system is a copper single crystal described by the embedded-atom method (EAM) 
potential. As shown in Fig.[T{a), the simulation cell is subjected to a pure shear stress along 
[112]. The dislocation to be nucleated lies on the (111) plane and has the Burgers vector 
of a Shockley partial^, b p = [1 12]/ 6. Fig. H^b) shows the shear stress-strain relationship of 
the perfect crystal at different temperatures (before dislocation nucleation). 

In this work, we predict the nucleation rate based on the Becker-Doring (BD) theory, 
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FIG. 2. (a) Free energy of dislocation loop during homog eneous nucleation at T — 300 K, cr X y — 
2.16 GPa (jxy = 0.135) from umbrella sampling, (b) Size fluctuation of critical nuclei from MD 
simulations. 



which expresses the nucleation rate per nucleation site as, 



rBD 



fc 1 eX P 



(2) 



k B T 

where / c + is the molecular attachment rate, and V is the Zeldovich factor (see Methods). 
The BD theory and TST only differs in the frequency prefactor. Whereas TST neglects 
multiple recrossing over the saddle point by a single transition trajectory^, the recrossing is 
accounted for in the BD theory through the Zeldovich factor. 

First, we establish the validity of the BD theory for dislocation nucleation by comparing 
it against direct MD simulations at a relatively high stress a = 2.16 GPa (7 = 0.135) at 
T = 300K, which predicts J MD = 2.5 x 10 8 s" 1 (see Methods). The key input to the BD 
theory is the activation Helmholtz free energy F C (T, 7), which is computed by umbrella 
sampling. The umbrella sampling is performed in Monte Carlo simulations using a bias 
potential as a function of the order parameter n, which is chosen as the number of atoms 
inside the dislocation loop (see Methods). 

Fig. |2](a) shows the free energy function F(n) obtained from umbrella sampling for the 
specified (T, 7) condition. The maximum of F(n) gives the activation free energy F c = 
0.53 ± 0.01 eV and the critical nucleus size n c = 36. The Zeldovich factor—, T = 0.055, is 
obtained from T = ( 27r k r ) where rj = — d 2 F(n)/dn 21 



in=n c 



Using the configurations collected from umbrella sampling with n = n c as initial condi- 



tions, MD simulations give the attachment rate /+ = 5.3 x 10 14 s _1 (see Methods). Because 
the entire crystal is subjected to uniform stress, the number of nucleation sites is the total 
number of atoms, iV atom = 14, 976. Combining these data, the Becker-Doring theory predicts 
the total nucleation rate to be A r atom / BD = 6.2 x 10 8 s -1 , which is within a factor of 3 of 
the MD prediction. The difference between the two is comparable to our error bar. This 
agreement is noteworthy because no adjustable parameters (such as the frequency prefactor) 
is involved in this comparison. It shows that the Becker-Doring theory and our numerical 
approach are suitable for the calculation of dislocation nucleation rate. 

We now examine the dislocation nucleation rate under a wide range of temperature and 
strain (stress) conditions relevant for experiments and beyond the limited time scale of MD 
simulations. Fig. [3]^a) shows the activation Helmholtz free energy F C (T, 7) as a function of 
7 at different T. The zero temperature data is obtained a minimum-energy-path (MEP) 
search using a modified version of the string method, similar to that used in^^i. The 
downward shift of F c curves with increasing T is the signature of the activation entropy 
S c - Fig. [3]^c) plots F c as a function of T at 7 = 0.092. For T < 400K, the data closely 
follow a straight line, whose slope gives ^(7) = 9 ks- This activation entropy contributes 
a significant multiplicative factor, exp(5 , c //cs) ~ 10 4 , to the absolute nucleation rate, and 
cannot be ignored. 

What causes this rapid drop of activation free energy with temperature? Thermal ex- 
pansion and vibrational entropy are two candidate mechanisms. To examine the effect of 
thermal expansion, we performed zero temperature MEP search at 7 = 0.092, but with other 
strain components fixed at the equilibrated values at T = 300 K. This approach is similar 
to the quasi-harmonic approximation (QHA)^>^ often used in free energy calculations in 
solids, except that, unlike QHA, the vibrational entropy is completely excluded here. The 
resulting activation energy, E c = 2.04 eV, is indistinguishable from the activation free energy 
F c = 2.05 ± 0.01 at T = 300 K computed from umbrella sampling. Because atoms do not 
vibrate in the MEP search, this result shows that the dominant mechanism for the large 
^(7) is thermal expansion, while the contribution from vibrational entropy is negligible. 
As temperature increases, thermal expansion pushes neighboring atoms further apart and 
weakens their mutual interaction. This expansion makes crystallographic planes easier to 
shear and significantly reduces the free energy barrier for dislocation nucleation. In the 
widely used harmonic approximation of TST, the activation entropy is often attributed to 
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FIG. 3. Activation free energy for homogeneous dislocation nucleation in copper, (a) F c as a 
function of shear strain 7 at different T. The data for T = 400K and T = 500K are shown in SI 
appendix, (b) G c as a function of shear stress a at different T. Squares represent umbrella sampling 
data and lines represent zero temperature MEP search results using simulation cells equilibrated at 
different temperatures, (c) F c as a function of T at 7 = 0.092. Circles represent umbrella sampling 
data and dashed line represent a polynomial fit. (d) G c as a function of T at a = 2.0 GPa from 
polynomial fit. 

the vibrational degrees of freedom as u exp(S c /k B ) = (YliLi U T) / (Yli^i 1 v t)i wnere V T anc ^ 
uf are the positive normal frequencies around the local energy minimum and activated state, 
respectively^^. However, here we see that 5 , c (7) arises entirely from the anharmonic effect 
for dislocation nucleation. At T = 400K and T = 500K, we observe significant differences 
between F c computed from umbrella sampling and E c computed from MEP search in the 
expanded cell. These difference must also be attributed to anharmonic effects. 

While it is easier to control strain 7 than stress a in atomistic simulations, it is usually 
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easier to apply stress in experiments, and experimental results are often expressed as a 
function of a and T. To bridge between simulations and experiments, it is important to 
establish a connection between the constant-stress and constant-strain ensembles. In the 
constant-strain ensemble, the system is described by the Helmholtz free energy F(n, T, 7) 
where n is the size of the dislocation loop and the activation Helmholtz free energy is defined 
as F C (T, 7) = F(n c ,T, 7) — F(n = 0,T, 7). In the constant-stress ensemble, the system is 
described by the Gibbs free energy G(n, T, a), from the Legendre transform G = F — a^V, 
with a = V~ 1 dF/d r y\ rit T- Similarly, G c (T,cr) = G{n c ,T,a) — G(n=0,T,a). We have proved 
that G C (T, a) = F C (T, 7) in the thermodynamic limit of V — > 00, when a and 7 satisfies the 
stress-strain relation of the perfect crystal, (7(7, T). The difference between F c and G c when 
cr = a(T, 7) is of the order OiV' 1 ). The details of the proof will be published separately. 

Combining the activation Helmholtz free energy F C (T, 7) shown in Fig.|3](a) and the stress- 
strain relations shown in Fig.[]Jb), we obtain the activation Gibbs free energy G C (T, a), which 
is shown in Fig. E^b). We immediately notice that the curves at different temperatures are 
more widely apart in G c (T,a) than that in F C (T, 7), indicating a much larger activation 
entropy in the constant-stress ensemble. For example, Fig. [3](d) plots G c as a function of T 
at <t = 2.0 GPa, from which we can obtain an averaged activation entropy of S c (a) = 48 k% 
in the temperature range of [0,300K]. This activation entropy contributes a multiplicative 
factor of exp(S' c (cr)/A; j B) ~ 10 20 to the absolute nucleation rate. 

The dramatic difference between S c (j) and S c (a) may seem surprising. Indeed, they are 
sometimes used interchangeably^^, although the conceptual difference between the two 
has been pointed out in the context of chemical reactions^2&. It is well known that the 
entropy is independent of the ensemble of choice, i.e. S(n, T, 7) = dF(n,T,-f)/dT\ n ^ and 
S(n,T,a) = dG(n,T,'j)/dT\ nja equal to each other as long as a = V^dF/d'yln^, which is 
true by definition. At the same time, the activation entropy is just the entropy difference 
between the activated state and the metastable state, i.e. S C (T, 7) = S(n c ,T, , ~f) — S(n= 
0,T, 7) and S C (T, a) = S(n c ,T,a) — S(n = 0,T,a). If the entropies in two ensembles can 
equal each other, it may seem puzzling how the activation entropies can be different. 

The resolution of this apparent paradox is that under the constant applied stress, the 
nucleation of a dislocation loop causes a strain increase, i.e. a(n=0, T, 7) = a(n c , T, 7 + ), with 
7 + > 7. Based on this result, one can show that the difference in the activation entropies, 
AS C = S c (a) — 5 C (7), equals S(n=0,T, 7 + ) — S(n=0,T, 7), which is the entropy difference 
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of the perfect crystal at two slightly different strains. We can further show that AS C = 
— Q c (a)da/dT\j, where Q c = —dG c / da\x is the activation volume and —da/dT\'-f describes 
the thermal softening effect. Hence AS C is always positive for nucleation processes in solids 
driven by shear stress. In the case of homogeneous dislocation nucleation AS C as large as 
39 ks is observed for homogeneous dislocation nucleation, which is mainly caused by its large 
activation volume Q c . The numerical results enable us to examine the approximation^ based 
on the so-called thermodynamic "compensation law"—, which states that the activation 
entropy is proportional to the activation enthalpy (or energy). We find that Sd'y) can be 
roughly approximated by E c {pj)/T* with T* « 3000 K while S c (a) is not proportional to 
H c {a) (see SI appendix). 

To assess the applicability of these conclusions in heterogeneous nucleation, we studied 
dislocation nucleation from the corner of a [001]-oriented copper nanorod with {100} side 
surfaces under axial compression (see Methods). Fig. IU(b) plots the activation free energy 
barrier as a function of axial compressive stress er, which shows significant reduction of the 
activation free energy with temperature. For example, at the compressive elastic strain of 
e = 0.03, the compressive stress is a = 1.50 GPa at T = K. The activation entropy S c (e) at 
this elastic strain equals 9&g, whereas the activation entropy S c (a) at this stress equals 17kg. 
Fig. H(c) plots the contour lines of the predicted dislocation nucleation rate (per nucleation 
site) as a function of T and a. To show the physical effect of the large activation entropies, 
the dashed lines plot the rate predictions if the effect of S c (a) were completely neglected. 
Significant deviations between the two sets of contour lines are observed, especially for 
T > 300 K and a < 1.5 GPa. For example, at T = 300 K and a = 1.5 GPa (where 
a thick and a thin contour line cross), the neglect of activation entropy would cause an 
underestimate of the nucleation rate by 10 orders of magnitude. 

In summary, we have shown that the Becker-Doring theory combined with free energy 
barriers determined by umbrella sampling can accurately predict the rate of homogeneous 
dislocation nucleation. In both homogeneous and heterogeneous dislocation nucleation, a 
large activation entropy at constant elastic strain is observed, and is attributed to the 
weakening of atomic bonds due to thermal expansion. An even larger activation entropy 
is observed at constant stress, due to thermal softening. Both effects are anharmonic in 
nature, and emphasize the need to go beyond harmonic approximation in the application of 
rate theories in solids. We believe our methods and the general conclusions are applicable 
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FIG. 4. (a) Heterogeneous dislocation nucleation in a copper nanorod under compression, visualized 
by Atomeye^. (b) G c as a function compressive stress a at different T. (c) Contour lines of 
dislocation nucleation rate per site / as a function of T and a. The predictions with and without 
accounting for activation entropy S c (a) are plotted in thick and thin lines, respectively. The 
nucleation rate of I ~ 10 6 s _1 per site is accessible in typical MD time scales while nucleation rate 
of I ~ 10~ 4 -10~ 9 is accessible in typical experimental time scales, depending on the number of 
nucleation sites. 

to a wide range of nucleation processes in solids that are driven by shear stress, including 
cross slip, twinning and martensitic phase transformation. 

Methods: 

Molecular Dynamics: The simulation cell for homogeneous dislocation nucleation has 
dimension 8 [112] x 6[111] x 3[110]. Periodic boundary conditions (PBC) are applied 
to all 3 directions. To reduce artifacts from periodic image interactions, the applied 
stress is always large enough so that the diameter of critical dislocation loop is smaller 
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than half the width of the simulation cell. 

The shear strain 7 is the x-y component of the engineering strain. The follow- 
ing procedure is used to obtain pure shear stress-strain curve shown in Fig. [H(b). 
At each temperature T and shear strain 7^, a series of 2 ps MD simulations un- 
der the canonical, constant temperature-constant volue (NVT) ensemble are per- 
formed. After each simulation, all strain components except j xy are adjusted ac- 
cording to the average Virial stress until a xy is the only non-zero stress component. 
The shear strain is then increased by 0.01 and the process repeats until the crystal col- 
lapses spontaneously. The shear stress-strain data are fitted to a polynomial function, 

To obtain average nucleation time at a xy = 2.16 GPa (7 = 0.135) at 300 K, we per- 
formed 192 independent MD simulations using the NVT ensemble with random initial 
velocities. Each simulation runs for 4 ns. If dislocation nucleation occurs during this 
period, the nucleation time is recorded. This information is used to construct the 
function P s (t), which is the fraction of MD simulation cells in which dislocation nucle- 
ation has not occurred at time t. P s (t) can be well fitted to the form of exp(— I MD t) 
to extract the nucleation rate J MD . 

To compute the attachment rate /+, we collect from umbrella sampling an ensemble 
of 500 atomic configurations for which n = n c , and run MD simulations using each 
configuration as an initial condition. The initial velocities are randomized according 
to Boltzmann's distribution. The mean square change of the loop size, (An 2 (t)), as 
shown in Fig. EJ^b), is fitted to a straight line, 2f+t, in order to extract f+— . 

Free energy barrier calculations: The reaction coordinate n is defined for each atomic 
configurations in the following way. An atom is labelled as "slipped" if its distance from 
any of its original nearest neighbors has changed by more than the critical distance 
d^.. We choose d c = 0.33A, 0.38A and 0.43A for T < 400 K, T = 500 K and 
T = 600 K, respectively. The "slipped" atoms are grouped into clusters; two atoms 
belong to the same cluster if their distance is less than cutoff distance r c (3.4A). The 
reaction coordinate n is the number of atoms in the largest cluster divided by two. 

To perform umbrella sampling, a bias potential fcgT (n — n) 2 is superimposed on the 
EAM potential, where T = 40 K and n is the center of the sampling window. We 
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choose T empirically so that the width of the sampling window on the n-axis is about 
10. The activation Helmholtz free energy for homogeneous nucleation data can be 
fitted very well by a polynomial function, F C (7,T) = ELo Ej=o ^j7 l ^ J lTy ^ ne range 
of 0.09 < 7 < 0.12 and < T < 500K. 

For heterogeneous dislocation nucleation, the size of the copper nanorocr^ is 15 [100] x 
15[010] x 20[001] with PBC along [001]. The activation Gibbs free energy for hetero- 
geneous nucleation is fitted to an empirical form G c (a,T) = A((a/a ) p — l) q — Ba l T . 

The error bar of activation free energies is about 0.5/cb^- Hence, the error bar of 
activation entropies is about 0.5/cb- 
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